home *** CD-ROM | disk | FTP | other *** search
/ HPAVC / HPAVC CD-ROM.iso / 3DGPL.ZIP / 3DGPL / TEXT / 3.TXT < prev    next >
Text File  |  1995-06-22  |  49KB  |  1,179 lines

  1.  
  2.  
  3.  
  4.  
  5.  
  6.                              2-D rasterization.
  7.                             --------------------
  8.  
  9.  
  10.  
  11.  
  12.  
  13.  
  14.                                 Re graphics: A picture is worth 10K words -
  15.                                 but only those to describe the picture.
  16.                                 Hardly any sets of 10K words can be adequately
  17.                                 described with pictures.
  18.  
  19.  
  20.  
  21.  
  22.  
  23.  
  24. Rasterization, the term itself refers to the technology we are using to
  25. represent images: turning it into mosaic of small square pixel cells set
  26. to different colours. The algorithms performing just that are covered in
  27. this section.
  28.  
  29. A dot.
  30. ------
  31. Rasterizing dot is very straightforward, I am not even sure it is
  32. warranted to call it rasterization at all. All it involves is just setting
  33. some location in the screen memory to particular colour. Let's recall
  34. that the very first location of the colourmap describes top rightmost
  35. pixel, that determines the system of references which is convenient in
  36. this case with X axis directed to the right and Y axis directed downward.
  37.  
  38.  
  39.                                HW_X_SIZE
  40.  
  41.                    (0,0)+-+-+-+- ... -+------->  X
  42.                         | | |         |
  43.                         +-+-+         +
  44.                         |
  45.             HW_Y_SIZE    ...
  46.  
  47.                         |
  48.                         +-+-      +-+
  49.                         | |       |*|<------y
  50.                         +-+-+     +-+
  51.                         |          ^
  52.                         |          |
  53.                         V          |
  54.                       Y            x
  55.  
  56.  
  57.                                 pic 3.1
  58.  
  59.  
  60. Assuming that each pixel is represented by one byte, and that the dimensions
  61. of the output surface is HW_SCREEN_X_SIZE by HW_SCREEN_Y_SIZE it is the
  62. y*HW_SCREEN_X_SIZE+x location we have to access in order to plot a dot
  63. at (x,y). Here's some code to do just that:
  64.  
  65.  
  66. unsigned char *G_buffer;                    /* colourmap- offscreen buffer */
  67.  
  68. void G_dot(int x,int y,unsigned char colour)
  69. {
  70.  G_buffer[y*HW_SCREEN_X_SIZE+x]=colour;
  71. }
  72.  
  73. A line.
  74. -------
  75. Line rasterization is a very important peace of code for the library. It
  76. would also be used for polygon rendering, and few other routines. What we
  77. have to do is set to some specified colour all dots on the line's path.
  78. We can easily calculate coordinates of all dots belonging to this line:
  79.  
  80.                                          * (xend,yend) ---
  81.                                        /                 |
  82.                                       * (x,y)            dy
  83.                                     / |
  84.                                   /   |y                 |
  85.                 (xstart,ystart) *-----+          ---------
  86.                                    x
  87.                                          |
  88.                                 |        |
  89.                                 |- -dx- -|
  90.  
  91.                                   pic 3.2
  92.  
  93. Remembering a little bit of high (high?) school geometry, it can be
  94. seen that if:
  95.                               dx = xend-xstart
  96.                               dy = yend-ystart
  97. then:
  98.                       x     dx               x * dy
  99.                      --- = ----    and  y = --------
  100.                       y     dy                 dx
  101.  
  102. Since we want to set all pixels along the path we just have to take
  103. all Xs in the range of [xstart, xend] and calculate respective Y. There
  104. is a bit of a catch, we would have only as much as xend-xstart
  105. calculated coordinates. But if the Y range was actually longer the path
  106. we have just found won't be continuous pic 3.3. In this case we should have
  107. been calculating coordinates taking values from the longer range
  108. [ystart, yend] and getting respective X, pic 3.4 as:
  109.  
  110.                                   y * dx
  111.                              x = --------
  112.                                     dy
  113.  
  114.               ....*                             ....*
  115.               .....                             ....*
  116.               ...*.                             ...*.
  117.               .....                             ...*.
  118.               ..*..  yrange                     ..*..  yrange
  119.               .....                             ..*..
  120.               .*...                             .*...
  121.               .....                             .*...
  122.               *....                             *....
  123.  
  124.               xrange                            xrange
  125.  
  126.              pic 3.3                           pic 3.4
  127.  
  128.  
  129. What we just did is quite logical and easily programmable algorithm,
  130. but... it is not that commonly used, the reason has to do with the way
  131. we did the calculation, it involved division, and that is by far the
  132. slowest operation any computer does. What we would do instead doesn't
  133. involve a single division. It is called a Bresenham's integer based
  134. algorithm (Bresenham, 1965) and here we find all points using few 
  135. logical operations and integer addition. The idea is to find a bigger 
  136. range among [xstart,xend] and [ystart,yend] go along points in it and 
  137. have a variable saying when it is time to advance in the other range.
  138. Let's consider what is happening at some stage in line rasterization,
  139. let's assume X is a longer range (dx>dy) and that with the line we are
  140. rendering xstart<xend  and  ystart<yend:
  141.  
  142.                                   |
  143.                               |       | H(x+1,y+1)
  144.                              -+---|---*-
  145.                               |       |h/
  146.                         |     |   |   * I(x+1,real_y)
  147.                               |      /|
  148.                         + - - - - + - - - - + -
  149.                               |   /   |
  150.                         |     | / |   |l    |
  151.                               /       |
  152.                              -*---|---*-L(x+1,y)
  153.                               |P(x,y) |
  154.  
  155.  
  156.                                pic 3.5
  157.  
  158.  
  159. Pic 3.5 shows the situation where we just rendered a pixel P (Previous) at
  160. coordinates (x,y) and now have to make a decision where to go to next, to
  161. pixel L(x+1,y) (Lower) or H(x+1,y+1) (Higher). points P,H,L actually represent
  162. middles of respective pixels. Since actual line would pass somewhere in the
  163. middle between those two points we must plot the one whose centre is closer
  164. to the intersection point I(x+1,real_y). this can be measured by comparing
  165. segments h and l, resulted from the intersection: remembering dependency
  166. used in the previous method we can find that:
  167.  
  168.                                      dy * (x+1)
  169.                            real_y = ------------
  170.                                          dx
  171. and:
  172.                                                       dy * (x+1)
  173.                  h = y+1-real_y     =>     h = y+1 - ------------
  174.                                                           dx
  175.  
  176.                                                 dy * (x+1)
  177.                  l = real_y-y       =>     l = ------------ - y
  178.                                                     dx
  179. Now, we are interested in comparing l and h:
  180.  
  181.                                      dy * (x+1)
  182.                           l - h = 2 ------------  - 2y-1
  183.                                         dx
  184.  
  185.  
  186. So, if l-h>0 it means that l>h the intersection point L is closer to point
  187. H and the later should be rendered, otherwise L should be rendered. let's
  188. multiply both sides by dx:
  189.  
  190.                        dx(l - h)= 2dyx + 2dy - 2dxy - dx
  191.  
  192. Now, since dx assumed bigger then 0, signs of (l-h) and dx(l-h) would be
  193. the same. Let's call dx(l-h) as d and find it's sign and value at some step
  194. i and next step i+1:
  195.  
  196.                       d    = 2dyx    -2dxy    + 2dy - dx
  197.                        i         i-1      i-1
  198.  
  199.                       d    = 2dyx    -2dxy    + 2dy - dx
  200.                        i+1       i        i
  201.  
  202. We can also find the initial value d, in the very first point since before
  203. that x==0 and y==0:
  204.  
  205.                                d  = 2dy-dx
  206.                                 1
  207.                               -------------
  208.  
  209. Since sign of d determines which point to move next (H or L) lets find
  210. value of d at some step i assuming we know it's value at i-1 from
  211. previous calculations:
  212.  
  213.               d   - d  = 2dyx - 2dxy - 2dyx    + 2dxy
  214.                i+1   i       i      i      i-1       i-1
  215.  
  216.               d   - d  = 2dy(x - x   ) - 2dx(y - y   )
  217.                i+1   i        i   i-1         i   i-1
  218.  
  219. Depending on the decision taken in the previous point value of d in the
  220. next one would be either:
  221.  
  222.                when H was chosen since x - x    =1 and y - y   =1
  223.                                          i   i-1         i   i-1
  224.  
  225.                                 d   - d = 2dy - 2dx
  226.                                  i+1   i
  227.  
  228.                                 d   = d + 2dy - 2dx
  229.                                  i+1   i
  230.                                ---------------------
  231.  
  232. or:
  233.  
  234.                when L was chosen since x - x   =1 and y - y   =0
  235.                                          i   i-1        i   i-1
  236.  
  237.                                 d   -d = 2dy
  238.                                  i+1  i
  239.                                 d   =d + 2dy
  240.                                  i+1  i
  241.                                --------------
  242.  
  243. This may seem a little complicated but the implementation code certainly
  244. doesn't look this way. We just have to find the initial value of d, and
  245. then in the loop depending on it's sign add to it either 2dy-2dx, or
  246. just 2dy. since those are constants, they can be calculated before the
  247. actual loop. In the proof above we have assumed xstart<yend ystrat<yend,
  248. however, in real life we can't guarantee that. The way we can take care
  249. of it, is just changing increments to decrements when above assumptions
  250. don't hold. We also have to remember that in the loop when L point is
  251. picked it is only value belonging to the bigger range which is incremented
  252. (decremented) whereas when picking H point both X and Y have to be changed
  253. since that is when we are advancing along both axes.
  254.  
  255.  
  256. void G_line(int x,int y,int x2,int y2,unsigned char colour)
  257. {
  258.  int dx,dy,long_d,short_d;
  259.  int d,add_dh,add_dl;
  260.  register int inc_xh,inc_yh,inc_xl,inc_yl;  
  261.  register int i;                            
  262.  
  263.  dx=x2-x; dy=y2-y;                          /* ranges */
  264.  
  265.  if(dx<0){dx=-dx; inc_xh=-1; inc_xl=-1;}    /* making sure dx and dy >0 */
  266.  else    {        inc_xh=1;  inc_xl=1; }    /* adjusting increments */
  267.  if(dy<0){dy=-dy; inc_yh=-1; inc_yl=-1;}
  268.  else    {        inc_yh=1;  inc_yl=1; }
  269.  
  270.  if(dx>dy){long_d=dx; short_d=dy; inc_yl=0;}/* long range,&making sure either */
  271.  else     {long_d=dy; short_d=dx; inc_xl=0;}/* x or y is changed in L case */
  272.  
  273.  d=2*short_d-long_d;                        /* initial value of d */
  274.  add_dl=2*short_d;                          /* d adjustment for H case */
  275.  add_dh=2*short_d-2*long_d;                 /* d adjustment for L case */
  276.  
  277.  for(i=0;i<=long_d;i++)                     /* for all points in longer range */
  278.  {
  279.   G_dot(x,y,colour);                        /* rendering */
  280.  
  281.   if(d>=0){x+=inc_xh; y+=inc_yh; d+=add_dh;}/* previous point was H type */
  282.   else    {x+=inc_xl; y+=inc_yl; d+=add_dl;}/* previous point was L type */
  283.  }
  284. }
  285.  
  286.  
  287. As you can see this algorithm doesn't involve divisions at all. It does
  288. have multiplication by two, but almost any modern compiler would optimize
  289. it to 1 bit left shift - a very cheap operation. Or if we don't have faith
  290. in the compiler we can do it ourselves. Lets optimize this code a little
  291. bit. Whenever trying to optimize something we first have to go after
  292. performance of the code in the loop since it is something being done
  293. over and over. In this source above we can see the function call to G_dot,
  294. let's remember what is inside there,.... oh, it is y*HW_SCREEN_X_SIZE+x a
  295. multiplication... an operation similar in length to division which
  296. we just spent so much effort to avoid. Lets think back to the
  297. representation of the colourmap, if we just rendered a point and
  298. now want to render another one next to it how their addresses in the
  299. colourmap would be different? pic 3.6
  300.  
  301.                             +-+-+-
  302.                            A|*->|B
  303.                             +|+-+-
  304.                            C|V| |
  305.                             +-+-
  306.                             | |
  307.  
  308.                             pic 3.6
  309.  
  310. Well, if it is along X axis that we want to advance we just have to add 1,
  311. to the address of pixel A to get to pixel B, if we want to advance along
  312. Y we have to add to the address of A number of bytes in the colourmap
  313. separating A and C. this number is exactly HW_SCREEN_X_SIZE that is length
  314. of one horizontal line of pixels in memory. Now, using that, let's try
  315. instead of coordinates (x,y) to calculate it's address in the colourmap:
  316.  
  317.  
  318. void G_line(int x,int y,int x2,int y2,unsigned char colour)
  319. {
  320.  int dx,dy,long_d,short_d;
  321.  int d,add_dh,add_dl;
  322.  int inc_xh,inc_yh,inc_xl,inc_yl;
  323.  register int inc_ah,inc_al;          
  324.  register int i;                      
  325.  register unsigned char *adr=G_buffer;
  326.  
  327.  dx=x2-x; dy=y2-y;                          /* ranges */
  328.  
  329.  if(dx<0){dx=-dx; inc_xh=-1; inc_xl=-1;}    /* making sure dx and dy >0 */
  330.  else    {        inc_xh=1;  inc_xl=1; }    /* adjusting increments */
  331.  if(dy<0){dy=-dy; inc_yh=-HW_SCREEN_X_SIZE;
  332.                   inc_yl=-HW_SCREEN_X_SIZE;}/* to get to the neighboring */
  333.  else    {        inc_yh= HW_SCREEN_X_SIZE; /* point along Y have to add */
  334.                   inc_yl= HW_SCREEN_X_SIZE;}/* or subtract this */
  335.  
  336.  if(dx>dy){long_d=dx; short_d=dy; inc_yl=0;}/* long range,&making sure either */
  337.  else     {long_d=dy; short_d=dx; inc_xl=0;}/* x or y is changed in L case */
  338.  
  339.  inc_ah=inc_xh+inc_yh;
  340.  inc_al=inc_xl+inc_yl;                      /* increments for point adress */
  341.  adr+=y*HW_SCREEN_X_SIZE+x;                 /* address of first point */
  342.  
  343.  d=2*short_d-long_d;                        /* initial value of d */
  344.  add_dl=2*short_d;                          /* d adjustment for H case */
  345.  add_dh=2*short_d-2*long_d;                 /* d adjustment for L case */
  346.  
  347.  for(i=0;i<=long_d;i++)                     /* for all points in longer range */
  348.  {
  349.   *adr=colour;                              /* rendering */
  350.  
  351.   if(d>=0){adr+=inc_ah; d+=add_dh;}         /* previous point was H type */
  352.   else    {adr+=inc_al; d+=add_dl;}         /* previous point was L type */
  353.  }
  354. }
  355.  
  356.  
  357. We can see from this code that although the complexity of the preliminary
  358. part of the algorithm went up a bit, the code inside the loop is much shorter
  359. and simpler. There is just one thing left to add in order to make it all
  360. completely usable - clipping that is making sure our drawing doesn't go
  361. outside the physical boundaries of the colourmap, but this is covered in
  362. the next chapter.
  363.  
  364.  
  365.  
  366. Ambient polygons
  367. ----------------
  368. What we have to do is to fill with some colour inside the polygon's edges.
  369. The easiest and perhaps fastest way to do it is by using "scan line"
  370. algorithm. The idea behind it is to convert a polygon into a collection
  371. of usually horizontal lines and then render them one at a time. The lines
  372. have to be horizontal only because the most common colourmap would have
  373. horizontal lines contingent in memory which automatically allows to render
  374. them faster then vertical lines since increment by 1 is usually faster
  375. then addition. There is one inherited limitation to the algorithm, because
  376. at each Y heights there would be just one horizontal line only concave
  377. polygons can be rendered correctly. It can be seen that non concave polygons
  378. might need two or more lines sometimes pic 3.7 (that's by the way, in case
  379. you wondered, determines the difference in definitions between concave and
  380. non-concave polygons):
  381.  
  382.                              *\
  383.                              |  \    /*
  384.                             y|---\  /-|
  385.                              |    \/  |
  386.                              *--------*
  387.  
  388.                                pic 3.7
  389.  
  390. How can we represent the collection of horizontal lines? obviously we
  391. just need to know at which Y heights it is, X coordinate of the point
  392. where it starts and another X coordinate of the point where it ends.
  393. So lets have one "start" and one "end" location for every Y possible
  394. on the screen. Since every line is limited by two points each belonging
  395. to certain edge, we can take one edge at a time find all points belonging
  396. to it and for each of them change the "start" and "end" at particular Y
  397. heights. So that at the end we would have coordinates of all scan lines
  398. in the polygon pic 3.8:
  399.  
  400.                  start end
  401.                 +-----+---+              1 2 3 4 5 6 7 8
  402.                 |  1  | 5 |- - - - - ->  *------\
  403.                 |  2  | 8 |               \------------*
  404.                 |  2  | 7 |- - - - - - ->  \----------/
  405.                 |  3  | 7 |                  \------/
  406.                 |  4  | 6 |- - - - - - - ->   *----*
  407.                 +-----+---+
  408.                                 pic 3.8
  409.  
  410. Initially we would put the biggest possible value in all locations of
  411. "start" array and the smallest possible in the locations of "end" array.
  412. (Those are by the way are defined in <limits.h> and it is not a bad
  413. practice to use them). Each time we've calculated the point on the edge
  414. we have to go to "start" and "end" locations at Y heights and if X of
  415. the point less then what we have in "start" write it there. Similar to
  416. that if this same X is bigger then value in "end" location we have to
  417. put it there too. Because of the way the arrays were initialized first
  418. point calculated for every Y height would change both "start" and "end"
  419. location. Let's look at the code to find Xs for each edge, so called
  420. scanning function:
  421.  
  422. void GI_scan(int *edges,int dimension,int skip)
  423. {
  424.  signed_32_bit cur_v[C_MAX_DIMENSIONS];     /* initial values for Z+ dims */
  425.  signed_32_bit inc_v[C_MAX_DIMENSIONS];     /* increment for Z+ dimensions */
  426.  signed_32_bit tmp;
  427.  int dx,dy,long_d,short_d;
  428.  register int d,add_dh,add_dl;              
  429.  register int inc_xh,inc_yh,inc_xl,inc_yl;  
  430.  int x,y,i,j;
  431.  int *v1,*v2;                               /* first and second verteces */
  432.  
  433.  v1=edges; edges+=skip; v2=edges;           /* length ints in each */
  434.  
  435.  if(C_line_y_clipping(&v1,&v2,dimension))   /* vertical clipping */
  436.  {
  437.   dx=*v2++; dy=*v2++;                       /* extracting 2-D coordinates */
  438.   x=*v1++; y=*v1++;                         /* v2/v1 point remaining dim-2 */
  439.   dimension-=2;
  440.  
  441.   if(y<G_miny) G_miny=y;
  442.   if(y>G_maxy) G_maxy=y;
  443.   if(dy<G_miny) G_miny=dy;
  444.   if(dy>G_maxy) G_maxy=dy;                  /* updating vertical size */
  445.  
  446.   dx-=x; dy-=y;                             /* ranges */
  447.  
  448.   if(dx<0){dx=-dx; inc_xh=-1; inc_xl=-1;}   /* making sure dx and dy >0 */
  449.   else    {        inc_xh=1;  inc_xl=1; }   /* adjusting increments */
  450.   if(dy<0){dy=-dy; inc_yh=-1; inc_yl=-1;}
  451.   else    {        inc_yh=1;  inc_yl=1; }
  452.  
  453.   if(dx>dy){long_d=dx;short_d=dy;inc_yl=0;} /* long range,&making sure either */
  454.   else     {long_d=dy;short_d=dx;inc_xl=0;} /* x or y is changed in L case */
  455.  
  456.   d=2*short_d-long_d;                       /* initial value of d */
  457.   add_dl=2*short_d;                         /* d adjustment for H case */
  458.   add_dh=2*short_d-2*long_d;                /* d adjustment for L case */
  459.  
  460.   for(i=0;i<dimension;i++)                  /* for all remaining dimensions */
  461.   {
  462.    tmp=v1[i]; tmp<<=G_P; cur_v[i]=tmp;      /* f.p. start value */
  463.    tmp=v2[i]-v1[i]; tmp<<=G_P;              /* f.p. increment */
  464.    if(long_d>0) inc_v[i]=tmp/long_d;        /* if needed (non 0 lines) */
  465.   }
  466.  
  467.   for(i=0;i<=long_d;i++)                    /* for all points in longer range */
  468.   {
  469.    if(x<G_start[y])                         /* further then rightmost */
  470.    {
  471.     G_start[y]=x;                           /* the begining of scan line */
  472.     for(j=0;j<dimension;j++)
  473.      G_rest_start[j][y]=cur_v[j];            /* all other dimensions */
  474.    }
  475.  
  476.    if(G_end[y]<x)                           /* further the leftmost */
  477.    {
  478.     G_end[y]=x;                             /* the end of scan line */
  479.     for(j=0;j<dimension;j++)
  480.      G_rest_end[j][y]=cur_v[j];              /* and for other dimension */
  481.    }
  482.  
  483.    if(d>=0){x+=inc_xh;y+=inc_yh;d+=add_dh;} /* previous dot was H type */
  484.    else    {x+=inc_xl;y+=inc_yl;d+=add_dl;} /* previous dot was L type */
  485.    for(j=0;j<dimension;j++)
  486.     cur_v[j]+=inc_v[j];                     /* for all other dimensions */
  487.   }
  488.  }
  489. }
  490.  
  491. As you can see this edge scanning function doesn't have much differences
  492. with our most basic line rendering code. The only thing we are updating the
  493. vertical size of polygon so that after all edges would go through this
  494. function we would have the correct vertical dimensions in G_miny and
  495. G_maxy. There is another difference: scan function is designed to process,
  496. other dimensions, that is it would interpolate X and similar to that
  497. other values specified in vertices, finding them for all Ys in range.
  498. We would need that when adding interpolative shading feature. As to the
  499. ambient polygon rendering, After scanning is over we would just have
  500. to render all horizontal scan lines in the range [G_miny,G_maxy]:
  501.  
  502. void G_ambient_polygon(int *edges,int length,unsigned char colour)
  503. {
  504.  int new_edges[G_MAX_TUPLES];               /* verticaly clipped polygon */
  505.  int new_length;                            /* although no edges there yet */
  506.  register unsigned char *l_adr,*adr;
  507.  register int beg,end,i;
  508.  
  509.  GI_boarder_array_init();                   /* initializing the arrays */
  510.  
  511.  new_length=C_polygon_x_clipping(edges,new_edges,2,length);
  512.  
  513.  for(i=0;i<new_length;i++)
  514.   GI_scan(&new_edges[i*2],2,2);             /* Searching polygon boarders */
  515.  
  516.  if(G_miny<=G_maxy)                         /* nothing to do? */
  517.  {
  518.   l_adr=G_buffer+G_miny*HW_SCREEN_X_SIZE;   /* addr of 1st byte of 1st line */
  519.  
  520.   for(; G_miny<=G_maxy; G_miny++,
  521.       l_adr+=HW_SCREEN_X_SIZE)              /* rendering all lines */
  522.   {
  523.    adr=l_adr+(beg=G_start[G_miny]);         /* addr of the current line */
  524.    end=G_end[G_miny];                       /* ends here */
  525.  
  526.    for(;beg<=end;beg++) *adr++=colour;      /* rendering single line */
  527.   }
  528.  }
  529. }
  530.  
  531. As to the optimization it doesn't appear we can do a lot without
  532. leaving comfortable boundaries of C. On the other hand when there's
  533. really nothing else we can do to speed things up, when all manuals
  534. are read and all friends are out of ideas let's do it, let's go
  535. for assembly. I don't think rewriting all functions this way is
  536. appealing. Moreover only rewriting the real execution bottleneck
  537. would give considerable speed increase. That's, perhaps, one of
  538. the reason why modern compilers have such a nice feature as inline
  539. assembly, allowing to add low-level code directly into C source.
  540.  
  541. Looking at the code above we can see this line rendering loop. This is
  542. an obvious candidate, since it executes fairly often and does quite a
  543. lot of iterations, especially compared with surrounding code.
  544.  
  545. Different C compilers have different provision for inline assembly.
  546. GNU C compiler has quite powerful syntax allowing to specify assembly
  547. code. It is:
  548.  
  549.      asm("assembly_command %1,%0":"=constraint" (result)
  550.                                  :" constraint" (argument),
  551.                                   " constraint" (argument),
  552.                                            ...
  553.                                  :" clobbered_reg", "clobbered_reg" ...
  554.         );
  555.  
  556. The assembly instruction is specified with aliases allowing
  557. to add names of C variables into it. it also allows to specify
  558. the CPU registers the command, we are inserting, might clobber.
  559. It is particularly important if we are expecting to compile with
  560. optimizing option or have register variables in the code. We
  561. don't want the program trying to use contents of a register we
  562. just destroyed by inline assembly instruction. If we are specifying
  563. the clobbered registers, on the other hand, compiler would make
  564. sure that no useful information is in this registers at the
  565. moment inline assembly code starts to execute.
  566.  
  567. WATCOM C compiler has different provision for inline assembly:
  568.  
  569. return_type c_function(parameter1,parametr2, ...);
  570. #pragma aux c_fucntion=        \
  571.            "assembly_command"  \
  572.            "assembly_command"  \
  573.            parm [reg] [reg]    \
  574.            value [reg]         \
  575.            modify [reg reg];
  576.  
  577. It is implemented by providing sort of a pseudo c-function which
  578. is treated by the compiler pretty much like a macro. "parm" specifies
  579. which registers take in parameters from function's argument list,
  580. "value"- from what register, if any, the return value is taken and
  581. "modify", specifies clobbered registers.
  582.  
  583. for DJGPP compiler code substituting line rendering loop:
  584.  
  585.    for(;beg<=end;beg++) *adr++=colour;      /* rendering single line */
  586.  
  587. may look like:
  588.  
  589.    asm("movl  %0,%%ecx" :: "g" (end):"%ecx");
  590.    asm("subl  %0,%%ecx" :: "g" (beg):"%ecx");
  591.    asm("movl  %0,%%edi" :: "g" (adr):"%edi");
  592.    asm("movl  %0,%%eax" :: "g" (colour):"%eax");
  593.    asm("cld");                              /* increment when stosb */
  594.    asm("rep");
  595.    asm("stosb %al,(%edi)");                 /* rendering line here */
  596.  
  597. for WATCOM:
  598.  
  599. void a_loop(unsigned char *adr,int end,int beg,int colour);
  600. #pragma aux a_loop=                      \
  601.         "sub ecx,esi"                    \
  602.         "cld"                            \
  603.         "rep stosb"                      \
  604.         parm [edi] [ecx] [esi] [al];
  605.  
  606. Not a terribly complicated code to write but something giving
  607. considerable speed increase. Then again we might have used a C
  608. library function memset for filling chunks of memory:
  609.  
  610.    void *memset(void *s,int c,size_t n);
  611.  
  612.    memset(l_adr,colour,end-beg);            /* filling scan line */
  613.  
  614. This function might have had just that code, we wrote, inside of
  615. it. let's not believe assumptions neither pro nor con and go check
  616. that out. with DJGPP or, in this matter, most UNIX compilers we can
  617. just extract and disassemble memset code from libc.a (standard
  618. c library) and see for ourselves. let's go:
  619.  
  620.      ar -x libc.a memset.o
  621.      objdump -d memset.o
  622.  
  623. We would see something like:
  624.  
  625. memset.o:     file format coff-go32
  626.  
  627. Disassembly of section .text:
  628. 00000000 <_memset> pushl  %ebp
  629. 00000001 <_memset+1> movl   %esp,%ebp
  630. 00000003 <_memset+3> pushl  %edi
  631. 00000004 <_memset+4> movl   0x8(%ebp),%edi
  632. 00000007 <_memset+7> movl   0xc(%ebp),%eax
  633. 0000000a <_memset+a> movl   0x10(%ebp),%ecx
  634. 0000000d <_memset+d> jcxz   00000011 <_memset+11>
  635. 0000000f <_memset+f> repz stosb %al,%es:(%edi)
  636. 00000011 <_memset+11> popl   %edi
  637. 00000012 <_memset+12> movl   0x8(%ebp),%eax
  638. 00000015 <_memset+15> leave
  639. 00000016 <_memset+16> ret
  640. 00000017 <_memset+17> Address 0x18 is out of bounds.
  641. Disassembly of section .data:
  642.  
  643. Very similar to what we wrote? Other compilers might have
  644. worse library code? it might be so, let's try WATCOM C,
  645. let's go:
  646.  
  647.      wlib clibs.lib *memset
  648.      wdisasm memset.obj
  649.  
  650. that's what we would see:
  651.  
  652. Module: memset
  653. Group: 'DGROUP' CONST,CONST2,_DATA,_BSS
  654.  
  655. Segment: '_TEXT' BYTE USE32  00000023 bytes
  656.  0000  57                memset          push    edi
  657.  0001  55                                push    ebp
  658.  0002  89 e5                             mov     ebp,esp
  659.  0004  8b 45 10                          mov     eax,+10H[ebp]
  660.  0007  8b 4d 14                          mov     ecx,+14H[ebp]
  661.  000a  8b 7d 0c                          mov     edi,+0cH[ebp]
  662.  000d  06                                push    es
  663.  000e  1e                                push    ds
  664.  000f  07                                pop     es
  665.  0010  57                                push    edi
  666.  0011  88 c4                             mov     ah,al
  667.  0013  d1 e9                             shr     ecx,1
  668.  0015  f2 66 ab                          repne   stosw
  669.  0018  11 c9                             adc     ecx,ecx
  670.  001a  f2 aa                             repne   stosb
  671.  001c  5f                                pop     edi
  672.  001d  07                                pop     es
  673.  001e  89 f8                             mov     eax,edi
  674.  0020  c9                                leave
  675.  0021  5f                                pop     edi
  676.  0022  c3                                ret
  677.  
  678. No disassembly errors
  679.  
  680. ------------------------------------------------------------
  681.  
  682. Is it not very similar to what we just did using inline assembly?
  683. WATCOM C code is even trying to be smarter then us by storing as
  684. much as it can by word store instruction, since every word store
  685. in this case is equivalent to two char stores, there should be
  686. twice less iterations. Of course there would be some time spent
  687. on the function call itself but most likely performance of the
  688. inline assembly code we wrote and memset call would be very very
  689. similar. It is just yet another example of how easy it might be
  690. in some cases to achieve results by very simple means instead of
  691. spending sleepless hours in front of the glowing box (unfortunately
  692. only sleepless hours teach us how to do things by simple means but
  693. that's the topic for another story only marginally dedicated to
  694. 3-D graphics).
  695.  
  696. By the way, in the provided source code there is an independent
  697. approach to fast memory moves and fills. That is functions to perform
  698. those are brought out into hardware interface as: HW_set_char ;
  699. HW_set_int; HW_copy_char; HW_copy_int. And their implementation
  700. may actually be either of the ways described above depending on
  701. the particular platform.
  702.  
  703.  
  704. Shaded polygons.
  705. ----------------
  706.  
  707. Shaded polygons are used to smooth appearance of faceted models,
  708. that is those where we approximate real, curved surfaces by a
  709. set of plane polygons. Interpolative or Gouraud shading (Gouraud, 1971), 
  710. we would discuss, is the easiest to implement, it is also not very 
  711. expensive in terms of speed. The idea is that we carry a colour 
  712. intensity value in every vertex of the polygon and linearly interpolate 
  713. those values to find colour for each pixel inside the polygon.
  714.  
  715. Finally this is the place where implemented N-dimensional scanning
  716. procedure comes in handy. Indeed colour intensity can be pretty
  717. much handled as just an extra dimension as far as edge scanning
  718. and clipping are concerned.
  719.  
  720.                              *I3
  721.                             /  \
  722.                            /      \
  723.                         I0*          \
  724.                            \           *I2
  725.                             \    I    /
  726.                           IT1*---*---*IT2
  727.                               \     /
  728.                                \   /
  729.                                 \ /
  730.                                  *I1
  731.  
  732.                               pic 3.9
  733.  
  734. We can obtain values on the left and right boarders of a polygon (IT1,
  735. IT2 on pic 3.9) by interpolating colour intensity value in every edge
  736. during edge scanning, just as "start/end' X values were found for every
  737. horizontal line. Afterwards when rendering certain horizontal scan line we
  738. can further interpolate "start" and "end" intensities for this line finding
  739. colour for pixels belonging to it (I on pic 3.9). Fixed point math is
  740. probably a very convenient way to implement this interpolation:
  741.  
  742. void G_shaded_polygon(int *edges,int length)
  743. {
  744.  int new_edges[G_MAX_SHADED_TUPLES];        /* verticaly clipped polygon */
  745.  int new_length;
  746.  register unsigned char *l_adr,*adr;
  747.  register int beg,end,i;
  748.  register signed_32_bit cur_c,inc_c;        /* current colour and it's inc */
  749.  
  750.  GI_boarder_array_init();                   /* initializing the array */
  751.  
  752.  new_length=C_polygon_x_clipping(edges,new_edges,3,length);
  753.  
  754.  for(i=0;i<new_length;i++)
  755.   GI_scan(&new_edges[i*3],3,3);             /* Searching polygon boarders */
  756.  
  757.  if(G_miny<=G_maxy)                         /* nothing to do? */
  758.  {
  759.   l_adr=G_buffer+G_miny*HW_SCREEN_X_SIZE;   /* addr of 1st byte of 1st line */
  760.  
  761.   for(; G_miny<=G_maxy; G_miny++,
  762.       l_adr+=HW_SCREEN_X_SIZE)              /* rendering all lines */
  763.   {
  764.    adr=l_adr+(beg=G_start[G_miny]);         /* addr of the current line */
  765.    end=G_end[G_miny];                       /* ends here */
  766.  
  767.    cur_c=G_0_start[G_miny];                 /* colour starts with this value */
  768.    inc_c=G_0_end[G_miny]-cur_c;
  769.    if(end>beg) inc_c/=end-beg+1;            /* f.p. increment */
  770.  
  771.    for(;beg<=end;beg++)                     /* render this line */
  772.    {
  773.     *adr++=cur_c>>G_P;                      /* rendering single point */
  774.     cur_c+=inc_c;                           /* incrementing colour */
  775.    }
  776.   }
  777.  }
  778. }
  779.  
  780. As you see code above looks not that much different from ambient polygon
  781. rendering source.
  782.  
  783. There is one small catch. I said that intensities can be handled pretty
  784. much as an extra dimension. In fact we can consider shaded polygon on
  785. the screen to take 2 space dimensions, and to have one colour dimension.
  786. But would all vertexes of this polygon belong to one plane in such 3-D space?
  787. If they would not, shading would look rather weird especially when this
  788. polygon would start to rotate. The common solution is limit polygons to just
  789. triangles. Why? well three points always belong to the same plane in
  790. 3-D space.
  791.  
  792.  
  793.  
  794. Textured polygons.
  795. ------------------
  796.  
  797. It might seem that we can extend interpolation technique used to
  798. calculate colour intensity for shading, onto texture rendering. Just
  799. keep texture X and Y in every vertex of the polygon, interpolate
  800. those two along edges and then along horizontal lines obtaining
  801. this way texture (X,Y) for every pixel to be rendered. This however 
  802. does not work... Or to be more exact it stops working when perspective 
  803. projection is used, the reason is simple: perspective transformation 
  804. is not linear.
  805.  
  806. One may ask, why did it work for interpolative shading, after all
  807. we were using perspective transformation all along? The answer is
  808. it did not work, but visual aspect of neglecting perspective effect
  809. for shading is quite suitable, neglecting it for texture rendering
  810. however is not. I believe it has to do with different image frequencies.
  811.  
  812. Just consider the following situation: a rectangular polygon is
  813. displayed on the screen so that it's one side is much closer
  814. to us then the other one, which is almost disappearing in the
  815. infinity: where do you think in the texture map lines A,B,C and D
  816. would be mapped from by linear interpolating method? how do you think
  817. the texture on the polygon would look like?
  818.  
  819.  
  820.             +\                       +--/---/+1
  821.            A|--\                    A|/  /   |
  822.            B|----\                  B|/      |
  823.            C|----/                  C|\    <------ Missed area
  824.            D|--/                    D|\  \   |
  825.             +/                       +--\---\+2
  826.  
  827.        Polygon on the                 Texture
  828.           Screen                        Map
  829.  
  830.  
  831.                          pic 3.10
  832.  
  833. Well it would look like if someone has carved triangular area
  834. marked "missed" on the picture above, and connected points 1 and
  835. 2 together, real ugly... In less dramatic cases texture rendered this
  836. way might look nicer, being perfect when effect of perspective
  837. transformation is negligible: all points are at almost the same
  838. distance from the viewer, or very far away from the viewing plane.
  839.  
  840. To get always nice looking results we would have to consider
  841. what is really happening when texture is being rendered:
  842.  
  843.  
  844.                                  u^                  Texture Space
  845.                                   | *T(u,v)
  846.                                   +---->
  847.                                  o     v
  848.  
  849.                          y ^
  850.                            |       z
  851.                            |  U \ /*W(x,y,z)/ V      World Space
  852.                            |    / \       /
  853.                            |  /     \   /
  854.                            |/         * O
  855.                            +--------------> x
  856.  
  857.                                 j^
  858.                                  |
  859.                                  |    *S(i,j)        Screen Space
  860.                                  |
  861.                                  +------->i
  862.  
  863.  
  864.                                    pic 3.11
  865.  
  866. Consider three spaces pictured above: the texture space, world space (the
  867. one before projection) and finally contents of the screen - screen or image
  868. space. We may think of them just as original polygon/texture map in case
  869. of texture space, polygon/texture map turned somehow in case of world
  870. space and finally same projected onto the screen.
  871.  
  872. If we know mapping of origin and two orthogonal vectors from texture
  873. space into the world space, we can express mapping of any point through
  874. them:
  875.  
  876.                            x = Ox + v*Vx + u*Ux
  877.                            y = Oy + v*Vy + u*Uy
  878.                            z = Oz + v*Vz + u*Uz
  879.  
  880. Where Vx...Uz are corresponding components of the respective vectors.
  881. Further point in the world space W(x,y,z) is being perspectively
  882. transformed into the screen space S(i,j):
  883.  
  884.                                  i = x / z
  885.                                  j = y / z
  886.  
  887. (the actual transformation used, would most likely also contain a
  888. multiplier - focus, but let's worry about particularities on some later
  889. stage). In order to perform mapping, on the other hand, we need u,v
  890. texture coordinates as function of screen coordinates i,j
  891.  
  892.                                  x = i * z
  893.                                  y = i * z
  894.  
  895. or:
  896.                   Ox + v*Vx + u*Ux = i * [ Oz + v*Vz + u*Uz ]
  897.                   Oy + v*Vy + u*Uy = j * [ Oz + v*Vz + u*Uz ]
  898.  
  899. trying to express u,v through i,j:
  900.  
  901.                     v*(Vx-i*Vz) + u*(Ux-i*Uz) = i*Oz - Ox
  902.                     v*(Vy-j*Vz) + u*(Uy-j*Uz) = j*Oz - Oy
  903.  
  904. further:
  905.                           (i*Oz - Ox) - u*(Ux - i*Uz)
  906.                      v = -----------------------------
  907.                                  (Vx - i*Vz)
  908.  
  909.                           (i*Oz - Ox) - v*(Vx - i*Vz)
  910.                      u = -----------------------------
  911.                                  (Ux - i*Uz)
  912. and:
  913.  
  914.      (i*Oz - Ox) - u*(Ux - i*Uz)
  915.     ----------------------------- * (Vy - j*Vz) + u*(Uy - j*Uz) = j*Oz - Oy
  916.               (Vx - i*Vz)
  917.  
  918.                      (i*Oz - Ox) - v*(Vx - i*Vz)
  919.     v*(Vy - j*Vz) + ----------------------------- * (Uy - j*Uz) = j*Oz - Oy
  920.                               (Ux - i*Uz)
  921.  
  922. or in nicer form:
  923.  
  924.           i*(Vz*Oy-Vy*Oz) + j*(Vx*Oz-Vz*Ox) + (Vy*Ox-Vx*Oy)
  925.       u = --------------------------------------------------
  926.           i*(Vy*Uz-Vz*Uy) + j*(Vz*Ux-Vx*Uz) + (Vx*Uy-Vy*Ux)
  927.  
  928.           i*(Uy*Oz-Uz*Oy) + j*(Uz*Ox-Ux*Oz) + (Ux*Oy-Uy*Ox)
  929.       v = --------------------------------------------------
  930.           i*(Vy*Uz-Vz*Uy) + j*(Vz*Ux-Vx*Uz) + (Vx*Uy-Vy*Ux)
  931.  
  932. At this point we have formulas of from where in the texture space the
  933. point in the screen space originates. There are several implementational
  934. problems, first few are simple, the actual perspective transformation
  935. is  i=x*focus/z rather then i=x/z, plus, we are performing screen centre
  936. translation so that screen space (0,0) is in the middle of the screen
  937. and not in the upper left corner. Hence original dependency should have
  938. been:
  939.  
  940.                i = (x * focus) / z + HW_SCREEN_X_CENTRE
  941.                j = (y * focus) / z + HW_SCREEN_Y_CENTRE
  942.  
  943.  
  944. And so finally, so called, reverse mapping formulas have to be amended 
  945. respectively:
  946.  
  947.                  ((i-HW_SCREEN_X_CENTRE)/focus) * (Vz* ...
  948.             u = -------------------------------------------
  949.                                 ...
  950.  
  951.                                 ...
  952.  
  953.  
  954. Another, a bit difficult part has to do with texture scaling.
  955. If the size of texture vectors was picked to be 1 that would
  956. assume no scaling, that is unit size along the polygon in the
  957. world space would correspond to unit size in the texture space.
  958. What we however want to do in a lot of instances is to map
  959. smaller texture onto a bigger polygon or the other way around,
  960. bigger texture onto a smaller polygon. Let's try scaling the
  961. texture space:
  962.  
  963.        T
  964.     *---------->
  965.    v| *        |         S
  966.     | | \          *--------------------->
  967.     | |    \   | v'| *                   |       v       T
  968.     | *-- ---*     | | \                       ----- = -----
  969.     V- - - - - |   | |   \               |       v'      S
  970.                    | |      \
  971.                    | |mapping \          |           v'*T
  972.                    | |of the      \              v = -----
  973.                    | |polygon      \     |             S
  974.                    | *---------------*
  975.                    V- - - - - - - - - - -|
  976.  
  977.                              pic 3.12
  978.  
  979.  
  980.  
  981. now we can put expression for v into direct mapping equations:
  982.  
  983.  
  984.                           Vx*T         Ux*T
  985.             x = Ox + v'* ------ + u'* ------
  986.                            S            S
  987.  
  988.                                ...
  989.  
  990. Vx and Ux multiplied by T in fact are just projections of vectors 
  991. enclosing the hole texture space, let's call them V'x == Vx*T, 
  992. U'x == Ux*T etc.
  993.  
  994.  
  995.                           V'x         U'x
  996.             x = Ox + v'* ----- + u'* -----
  997.                            S           S
  998.                                ...
  999.  
  1000.  
  1001. considering this reverse mapping equations would look like:
  1002.  
  1003.  
  1004.                 ((i-HW_SCREEN_X_CENTRE)/focus) * (V'z* ...
  1005.        u'= S * --------------------------------------------
  1006.                                 ...
  1007.  
  1008.                                 ...
  1009.  
  1010. This would allow us to scale texture on the polygon by changing S
  1011. multiplier.
  1012.  
  1013. Another problem has to do with point O - mapping of the origin of the
  1014. texture space into the world space. The thing is, if the mapping of this
  1015. point gets onto Z==0 plane mapping equations above no longer hold. Just 
  1016. due to the fact that perspective transformation is having singularity 
  1017. there. In general we deal with this problem of perspective transformation
  1018. by clipping the polygon against some plane in front of Z==0. And do we
  1019. really need mapping of exactly O of the texture space? Not necessarily,
  1020. it may be any point in texture space assuming we change a bit direct
  1021. mapping equations:
  1022.  
  1023.  
  1024.                       x' = Ox + (v'-v'0)*V'x + (u-u0)*U'x
  1025.                                    ...
  1026.  
  1027.  
  1028. where (v'0,u'0) are coordinates in the texture space of the point we
  1029. have a mapping for in the world space that is (Ox,Oy,Oz). And the
  1030. reverse mapping equations would be then:
  1031.  
  1032.                    ((i-HW_SCREEN_X_CENTRE)/focus) * (V'z* ...
  1033.            u'= S * ------------------------------------------- + u'0
  1034.                                    ...
  1035.  
  1036.                                    ...
  1037.  
  1038. How can we obtain a mapping of, now, any point from the texture
  1039. space into the world space? If we would associate with every
  1040. polygon's vertex besides just x,y,z also u,v of where this vertex
  1041. is in the texture, we can treat that as 5-D polygon and perform
  1042. clipping on it, obtaining vertexes with coordinates suitable
  1043. for perspective transformation, hence any of them would be fine 
  1044. for the mapping equations. (How to do clipping on polygons is 
  1045. covered in the next chapter).
  1046.  
  1047. Last thing, in order to render regular patterns, those repeating
  1048. themselves, we may want to introduce texture size parameter,
  1049. and make sure that if we want to access texture outside this
  1050. size, this access would wrap around to pick a colour from somewhere
  1051. inside the texture. This can be easily done if texture size
  1052. is chosen to be some power of 2, for in this case we can
  1053. perform range checking with just logical "and" operation.
  1054.  
  1055. However, using above equations for texture mapping appear to be
  1056. quite expensive, there are few divisions involved per each rendered
  1057. pixel. How to optimize that? The, perhaps, easiest way is: horizontal
  1058. line subdivision, that is calculating real texture mapping every N pixels
  1059. and linearly interpolating in between, this method is used in the
  1060. implementation below:
  1061.  
  1062.  
  1063. void G_textured_polygon(int *edges,int length,int *O,int *u,int *v,
  1064.             unsigned char *texture,int log_texture_size,
  1065.                            int log_texture_scale
  1066.                )
  1067. {
  1068.  int new_edges[G_MAX_SHADED_TUPLES];        /* verticaly clipped polygon */
  1069.  int new_length;
  1070.  signed_32_bit Vx,Vy,Vz;
  1071.  signed_32_bit Ux,Uy,Uz;                    /* extracting vectors */
  1072.  signed_32_bit x0,y0,z0;
  1073.  signed_32_bit ui,uj,uc;
  1074.  signed_32_bit vi,vj,vc;
  1075.  signed_32_bit zi,zj,zc;                    /* back to texture coeficients */
  1076.  signed_32_bit v0,u0;
  1077.  signed_32_bit xx,yy,zz,zzz;
  1078.  int xstart,xend;
  1079.  signed_32_bit txstart,tystart;
  1080.  signed_32_bit txend,tyend;
  1081.  signed_32_bit txinc,tyinc;
  1082.  register unsigned char *g_adr,*adr;
  1083.  register int i,x,y;
  1084.  signed_32_bit txtrmasc=(0x1<<(log_texture_size+G_P))-0x1;
  1085.  log_texture_scale+=G_P;
  1086.  
  1087.  GI_boarder_array_init();                   /* initializing the array */
  1088.  
  1089.  new_length=C_polygon_x_clipping(edges,new_edges,4,length);
  1090.  
  1091.  for(i=0;i<new_length;i++)
  1092.   GI_scan(&new_edges[i*4],2,4);             /* Searching polygon boarders */
  1093.  
  1094.  if(G_miny<=G_maxy)                         /* nothing to do? */
  1095.  {
  1096.   x0=O[0]; y0=O[1]; z0=O[2];
  1097.   u0=O[3]<<G_P; v0=O[4]<<G_P;               /* world point <-> texture point */
  1098.  
  1099.   Vx=v[0]; Vy=v[1]; Vz=v[2];
  1100.   Ux=u[0]; Uy=u[1]; Uz=u[2];                /* extracting vectors */
  1101.  
  1102.   ui=(Vz*y0)-(Vy*z0);
  1103.   uj=(Vx*z0)-(Vz*x0);
  1104.   uc=(Vy*x0)-(Vx*y0);
  1105.   vi=(Uy*z0)-(Uz*y0);
  1106.   vj=(Uz*x0)-(Ux*z0);
  1107.   vc=(Ux*y0)-(Uy*x0);
  1108.   zi=(Vy*Uz)-(Vz*Uy);
  1109.   zj=(Vz*Ux)-(Vx*Uz);
  1110.   zc=(Vx*Uy)-(Vy*Ux);                       /* back to texture */
  1111.  
  1112.   g_adr=G_buffer+G_miny*HW_SCREEN_X_SIZE;   /* addr of 1st byte of 1st line */
  1113.  
  1114.   for(; G_miny<=G_maxy; G_miny++,
  1115.       g_adr+=HW_SCREEN_X_SIZE)              /* rendering all lines */
  1116.   {
  1117.    xstart=G_start[G_miny];
  1118.    adr=g_adr+xstart;
  1119.    xstart-=HW_SCREEN_X_CENTRE;
  1120.    x=xstart;
  1121.    xend=G_end[G_miny]-HW_SCREEN_X_CENTRE;
  1122.    y=G_miny-HW_SCREEN_Y_CENTRE;
  1123.  
  1124.    xx=((y*uj)>>T_LOG_FOCUS)+uc;
  1125.    yy=((y*vj)>>T_LOG_FOCUS)+vc;
  1126.    zz=((y*zj)>>T_LOG_FOCUS)+zc;             /* valid for the hole run */
  1127.  
  1128.    if(( zzz=zz+((x*zi)>>T_LOG_FOCUS) )!=0)
  1129.    {
  1130.     txend=( (xx+ ((x*ui)>>T_LOG_FOCUS)) <<log_texture_scale )/zzz +u0;
  1131.     tyend=( (yy+ ((x*vi)>>T_LOG_FOCUS)) <<log_texture_scale )/zzz +v0;
  1132.    }
  1133.  
  1134.    for(;xstart<xend;)
  1135.    {
  1136.     x+=G_LINEAR; if(x>xend) x=xend;         /* size of linear run */
  1137.     txstart=txend;
  1138.     tystart=tyend;
  1139.  
  1140.     if(( zzz=zz+((x*zi)>>T_LOG_FOCUS) )!=0)
  1141.     {
  1142.      txend=( (xx+ ((x*ui)>>T_LOG_FOCUS)) <<log_texture_scale )/zzz +u0;
  1143.      tyend=( (yy+ ((x*vi)>>T_LOG_FOCUS)) <<log_texture_scale )/zzz +v0;
  1144.     }
  1145.  
  1146.     length=x-xstart;                        /* ends here */
  1147.     if(length!=0) { txinc=(txend-txstart)/length;
  1148.             tyinc=(tyend-tystart)/length;
  1149.           }
  1150.  
  1151.     for(;xstart<x;xstart++)                 /* linear run */
  1152.     {
  1153.      txstart&=txtrmasc; tystart&=txtrmasc;  /* wrap around */
  1154.      *adr++=*(texture+((tystart>>G_P)<<log_texture_size)+(txstart>>G_P));
  1155.      txstart+=txinc; tystart+=tyinc;        /* new position in the texture */
  1156.     }
  1157.    }
  1158.   }
  1159.  }
  1160. }
  1161.  
  1162.  
  1163.  
  1164.  
  1165. Other possible speed-up techniques are: area subdivision involving
  1166. 2-D interpolation, faking texture mapping with polynomials (later
  1167. has not very much to do with the true mapping equations described
  1168. here, however) and rendering texture along constant z lines (and
  1169. not along horizontal line) the advantage gained by the former is
  1170. that along every constant Z line perspective transformation is
  1171. linear ( perspective transformation is i=x/z if z==const we
  1172. can write it as i=coef*x where coef=1/z which is linear of course)
  1173.  
  1174. The function above might be extended with interpolative shading
  1175. as well. To do that special consideration has to be given to
  1176. palette's layout, that is where and how to keep different intensities
  1177. of the same colour. But once decided coding is quite straightforward.
  1178.  
  1179.                                    * * *